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Abstract. We present a method based on Dickson's lemma to compute the 
"approximate radical" of a zero dimensional ideal I in C[a;i, . . . ,Xm\ which 
has zero clusters: the approximate radical ideal has exactly one root in each 
cluster for sufficiently small clusters. Our method is "global" in the sense that 
it does not require any local approximation of the zero clusters: it reduces the 
problem to the computation of the numerical nuUspace of the so called "ma- 
trix of traces", a matrix computable from the generating polynomials of 7. 
To compute the numerical nuUspace of the matrix of traces we propose to use 
Gauss elimination with pivoting or singular value decomposition. We prove 
that if I has k distinct zero clusters each of radius at most e in the oo-norm, 
then k steps of Gauss elimination on the matrix of traces yields a submatrix 
with all entries asymptotically equal to e^. We also show that the (fc + l)-th 
singular value of the matrix of traces is proportional to . The resulting ap- 
proximate radical has one root in each cluster with coordinates which are the 
arithmetic mean of the cluster, up to an error term asymptotically equal to 
£^ . In the univariate case our method gives an alternative to known approx- 
imate square-free factorization algorithms which is simpler and its accuracy 
is better understood. 
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Introduction 

Let / C C[x] be a polynomial ideal in m variables x = [xi, . . . ,Xm] with roots 
Zi, . . . , Zfc G C™ of multiplicities ni, . . . , n^, respectively, and let / e C[x] be an 
ideal with clusters Ci , . . . , such that each cluster Ci has rii roots around Zi 
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within radius e in the oo-norm for i ~ 1, . . . , fc. We present an algorithm which 

computes an approximate radical of /, denoted by \/7, which has exactly one root 
for each cluster, and we show that such root corresponds to the arithmetic mean 
of the cluster. 

The method we present in the paper is "global" in the sense that we do 
not use any local information about the roots in the clusters, only the coeffi- 
cients of the system of polynomials defining /, and we return another system of 
polynomials where all near multiplicities are eliminated. In the univariate case 
such global algorithms are used for example in approximate factoring (see [25]). 
where the input polynomial needs to be "square-free" in the approximate sense. 
Previous global methods which handle univariate polynomials with clusters use 
approximate gcd computation and approximate polynomial division in order to ei- 
ther factor out the near multiplicities or to compute the approximate multiplicity 
structure and find the roots of the nearest polynomial with the given multiplicity 
structure [HI HU [551 [S] ■ The method we propose here offers an alternative algo- 
rithm to factor out near multiplicities, which is simpler, and the relation between 
the accuracy of the output and the size of the clusters is better understood. We 
describe separately our method applied to the univariate case, and illustrate its 
simplicity and accuracy. 

Our method is based on Dickson's lemma, which gives the Jacobson radical 
of a finite dimensional associative algebra over a field of characteristic via the 
vanishing of traces of elements in the algebra. An immediate application of Dick- 
son's lemma to the algebra C[x]// finds a basis for Vl/I by finding the nuUspace 
of the matrix of traces R, a matrix computable from the generating polynomials 
of / using either multiplication matrices or other trace computation methods, as 
described below. 

The main focus of the paper is to adapt the method based on Dickson's 
lemma to the case when the ideal / has clusters of roots. In the paper we assume 
that both C[x]// and C[x]// are finite dimensional over C and have the same 
basis B C C[x]. Note that if / is generated by a well-constrained system, then 
"almost all" perturbations J of / will satisfy our assumption, however our results 
are not limited to well-constrained systems only. On the other hand, the results we 
prove in this paper measure the accuracy of the output in terms of the size of the 
clusters, as opposed to the size of the perturbation of the generating polynomials 
of the ideal /. The extension of our method to handle perturbations which change 
the structure of the factor algebra and to understand the accuracy of the output 
in terms of the size of the coefficient perturbation is the subject of future research. 
The results in this paper can be summarized as follows: 

Given the basis B and the matrix of traces R associated to / and B, us- 
ing Gaussian elimination with complete pivoting (GECP) we give asymptotic es- 
timates of order for the "almost vanishing" entries in Uk, the partially row 
reduced matrix of R, as well as upper bounds for the coefficients of , where e 
is the radius of the clusters in the oo-norm. These bounds can be used to give a 
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threshold to decide on the numerical rank of R, and to indicate the relationship 
between the numerical rank and the size of the clusters. 

Alternatively, we show how our results for the GECP of the matrix of traces 
R imply asymptotic bounds on the singular values of R. We also obtain in this 
case that the "almost vanishing" singular values are proportional to the square of 
the size of the clusters. This implies that for the numerical rank determination of 
i?, computing its SVD works similarly as using GECP. 

Using a basis of the numerical nullspace of R (or possibly an extended ver- 
sion of it), we define a set of generating polynomials for the approximate radical 

ideal VT, or similarly, define a system of multiplication matrices M^^ , • . • , M'^^^ 
of C[x]/\/7 with respect to a basis B' . We prove that modulo the generat- 
ing polynomials of ^/l are consistent and have roots with coordinates which are 
the arithmetic means of the coordinates of the roots in the clusters, which implies 
that the matrices M^^ , . • ■ , M'^^ commute and their eigenvalues are the arithmetic 
means of the coordinates of the roots in the clusters, all modulo . In other words, 
our algorithm finds the coefficients of a polynomial system with roots which are 
the means of the clusters up to a precision of about twice as many digits as the 
radius of the clusters, assuming that the clusters are sufficiently small. 

Let us briefly mention some of the possible methods to compute the matrix 
of traces i?, although in the paper we do not elaborate on this aspect. As we 
shall demonstrate in the paper, the matrix of traces R is readily computable from 
a system of multiplication matrices of C[x]//, for example from Mx-^^, . . . , M^^, 
where Mx^ denotes the matrix of the multiplication map by Xi in C[x]// written 
in terms of the basis B. One can compute Mx^ using Grobner bases (see for example 
[9]), resultant and subresultant matrices [30l[8l[49j, Lazard's algorithm [221III], or 
by methods that combine these [35]. Thus, our algorithm reduces the problem of 
finding the eigenvalues of matrices Mx^ , . . . , Mx^ which have clustered eigenvalues 
to finding eigenvalues of the smaller matrices M'x^ , ■ ■ ■ , M'x^ with well separated 
eigenvalues. 

In certain cases, the matrix of traces can be computed directly from the 
generating polynomials of /, without using multiplication matrices. We refer to 
the papers [H [T71 [H [71 [5] for the computation of traces using residues and Newton 
sums, or [13] using resultants. 

Also, fast computation techniques like the "baby steps-giant steps" method 
[24] 146] [45] can be implemented to speed up the computation of all entries 
of the matrix of traces. As we prove in the paper, the entries of the matrix of 
traces R are continuous in the size e of the root perturbation around e = 0, unlike 
the entries of multiplication matrices which may have many accumulation points 
as e approaches zero. Therefore, avoiding the computation of the multiplication 
matrices has the advantage of staying away from the possible large computational 
errors caused by the discontinuity of their entries. 
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111 the multivariate case, most of the methods handhng clusters of roots in the 
literature are "local" in that they assume sufficiently close approximations for the 
clusters in question. Our algorithm, viewed as having the multiplication matrices 
as input, is closest to the approach in [301 llOj in that these papers also reduce 
the problem to the computation of the eigenvalues of a system of approximate 
multiplication matrices. Both of these papers propose to reorder the eigenvalues 
of the multiplication matrices to group the clusters together. For the reordering 
of the eigenvalues these papers compute approximations of the eigenvalues by 
either using the approach in [2] or using the univariate method of |21j . In contrast, 
our method reorders the eigenvalues of all multiplication matrices simultaneously 
without approximating the eigenvalues, grouping one eigenvalue from each of the 
clusters together in a way which facilitates the computation of the means of the 
clusters and the elimination of the rest of the nearly repeated eigenvalues. Another 
local method to handle near multiple roots is the "deflation" algorithm, studied 
in the works [HI [301 [551 [1S| , to replace the original system which had a near 
multiple root with another one which has the same root with multiplicity one, 
using an approximation of the root in question. Related to the deflation algorithm, 
in [471 mi E] methods are proposed to compute the multiplicity structure of 
a root locally in terms of the so called dual basis, and then computing good 
approximations for the individual roots in the cluster, assuming that either a 
near system with multiple roots is known, or a sufficient approximation of the 
multiple root is given. Additionally, methods for computing singular solutions 
of both polynomials and analytic functions using homotopy continuation can be 
found in [31 [Ml [35]. 

We also include here reference to some of the related methods for solving 
systems of polynomial equations with exact multiplicities: involving the computa- 
tion of dual bases |32l [31] [48] , or in the univariate (or bivariate) case, using Gauss 
maps |26j , or analyzing the structure of the multiplication matrices by transform- 
ing them to an upper triangular form |50 1 136| , [37]. Previous work using Dickson's 
Lemma to compute radical ideals in the exact case includes [HIS]. Also, [43] uses 
trace matrices in order to find separating linear forms deterministically. 

The present paper is the extended and unabridged version of the paper that 
appeared in [22]. 

Acknowledgements: 
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1. Preliminaries 

Let A be an associative algebra over a field F of characteristic 0. (See definition 
and basic properties of associative algebras in [191 142].) 

An element x € A is nilpotent if x"^ = for some positive integer m. 

An element x € A is properly nilpotent if xy is nilpotent for every y € A. 
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The radical of A, denoted Rad{A), is the set of properly nilpotent elements 
of A. The radical Rad{A) is an ideal of A. In commutative algebras nilpotent 
elements are properly nilpotent, hence for a commutative A the radical Rad{A) is 
simply the set of nilpotent elements in A. 

Throughout the paper we assume that A is finite dimensional over F. Fix a 
basis B = [bi, . . . ,bn] oi A (note that later we will need to fix the order of the 
elements in B, that is why we use vector notation). We call the multiplication 
matrix of a; £ A the matrix of the multiplication map 

rrix '■ A — > A 

[9] '-^ [xg] 

written in the basis B. It is easy to verify (cf. Page 8 in that the map x 1-^ 
is an algebra homomorphism, called regular representation from A to Mn{F). 

The trace of x, denoted Tr{x), is the trace of the matrix M^. It is independent 
of the choice of the basis. 

2. Matrix Traces and the Radical 

Our main construction is based on the following results describing the elements of 
the radical of an associative algebra A using traces of elements: 

Theorem 2.1 (Dickson [l^ pp. 106-107). An element x of an associative algebra A 
over a field F of characteristic is properly nilpotent if and only if Tr{xy) ~ 
for every y Cz A. 

Corollary 2.2 (Friedl and Ronyai (TH] p. 156). Let F be a field of characteristic 
and A a matrix algebra over F. Let B = [61, ... , 6„] be a linear basis of A over the 
field F. Then x G Rad[A) if and only if Tr{xbi) = 0, i = 1, ... ,71. 

We apply the above results to the special case of commutative algebras which 
are quotients of polynomial rings. Consider the system of polynomial equations 

f (x) = 

where f = {/i ,...,//} and each /; is a polynomial in the variables x = [xi , . . . , Xm] ■ 
Assume that the polynomials /i , . . . , /; have finitely many roots in C™ , which 
implies that the algebra A — C[x]// is finite dimensional, where / is the ideal 
generated by the polynomials in f . Denote the dimension of A over C by n and 
let B = [61, . . . , 6„] be a basis of A. By slight abuse of notation we denote the 
elements of the basis B which are in A and some fixed preimages of them in 
C[a;i, . . . , Xm] both by 5i, . . . , 6„. Let {zi, . . . , z„} C C™ be the set of common 
roots (not necessarily all distinct) of the polynomials in f . Using the multiplication 
matrices Mf associated to the elements f E A and the fact that 

RadiA) = Vl/I C C[x]// = A, 
we can reword Corollarv l2.2l in the following way: 
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Corollary 2.3. Letp e C[x] andp be the image ofp in A. Using the above notation, 
the following statements are equivalent: 

{i) peVi 
(a) p e Rad{A) 

{Hi) Tr{Mphj ) ~ for all j = 1, . . . , n. 

We can now use the previous corollary to characterize the radical of A as the 
nuUspace of a matrix defined as follows: 

Definition 2.4. The matrix of traces is the n x n symmetric matrix: 

where Mi,-bj is the multiplication matrix of bib j as an element in A in terms of the 
basis B = [&i , . . . , 6„] and Tr indicates the trace of a matrix. 

Corollary 2.5. An element 

n 

= ^ Ckbk 

k=l 

of the quotient ring A with basis i? = [foi , . . . , 6„] is in the radical of A if and only 
if [ci, . . . , Cn] is in the nullspace of the matrix of traces R. 

Proof. Corollarv l2. 31 states that an element r ~ X]fc=i ^fc^fe S A belongs to Rad{A) 
if and only if Tr{Mrb ) — 0, for all j — 1, . . . , n. From the linearity of both the 
multiplication map (see Proposition (4.2) in Chapter 2 of [121) and the traces of 
matrices we have that 

n 

Tr{Mrb,) = Y.CkTr{Nh,b,) 

k=l 

= [Ci, . . .,Cn]R[j] 

where R[j] is the j*'* column of the matrix of traces R. Therefore. Tr{Mrb ) — 
for all j = 1 , . . . , ft is 6(^uiv&,leiit to [ci, . ■ . , Cn 

□ 

Remark 2.6. Methods in the literature for computing the matrix of traces R are 
mentioned in the Introduction. One way to compute it is from the multiplication 
matrices Mb^bj ■ Note that in order to compute the matrices Mb^bj , j = ^, ■ ■ ■ ,n, 
it is sufficient to have M^^, k = 1, . . . , m, since if /i € C[xi, . . . , Xm] is a preimage 
of bibj € A, then we have 

= h{M^,,...,M^J. 



This is because the regular representation is a homomorphism of C-algebras, see 
also Corollary (4.3) in Chapter 2 of [12]. 
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Example 2.7. We consider the polynomial system fi = /2 = /a = 0, with 

/i —3^1 + 4:XiX2 — 6x1 + — 18x2 + 13 

/2 —x^ + 16x^X2 — 7x\ + 1182:ia;2 — 286x12:2 

+ 147x1 - X2 + 6x2 + X2 + 5 
/a — x'l + 10x1X2 — 5xi + 72x1X2 — 176x1X2 

+ 91x1 — X2 + 4x2 + X2 + 3 

These polynomials have two common roots: [1, 1] of multiplicity 3 and [—1,2] of 
multiplicity 2. 

We compute the multiplication matrices Mx-^ and with respect to the 

basis B — [1, xi, a;2, xia;2, x^], which are respectively 







1 

-1 2 

^ 3 

4 4 



(1) 



Here we used Chardin's subresultant construction to compute the multiplica- 
tion matrices. (See and |49] .) 

We now compute the matrix R using Definition \2.4\ and Remark[ 



R = 



5 1 

1 5 

7 -1 

-1 7 

5 1 



7 -1 

-1 7 

11 -5 

-5 11 

7 -1 



(2) 



The nullspace of R is generated by the vectors 

[1, -3, 0, 2, 0], [0, -4, 1, 3, 0]. [0, -3, 0, 2, 1]. 

By Corollarv \2.5\ we have that the radical of I = (/i,/2,/3) modulo I is 

^/l 1 1 — ^ ~ 3x1 -!- 2xiX2, —4x1 -f X2 + 3xiX2, —3x1 + 2xiX2 + x^^ . 

Note that the polynomials on the right hand side are in \fl . 

Assume that rank R ~ k. Once we know the n~k generators {rk+i, • ■ • , T'n} of 
the radical, we can obtain the multiplication matrices of the elements oi A/ Rad{A) — 
C [x] /\/7 by performing a change of basis on the multiplication matrices M^-^ , • ■ ■ , 
to the basis {ri, . . . ,rk, rk+i, . . . , r„} of A, where ri, . . . , r/j can be chosen arbi- 
trarily as long as {ri, . . . , r^, rk+i, ■ ■ ■ , ?■„} is linearly independent. Let M^^ be the 
multiplication matrix of the coordinate Xs in the basis [ri, . . . , r„]. Then the k x k 
principal submatrix 

is the multiplication matrix of Xs in A/Rad{A) — C[x]/'\/7 with respect to the 
basis [ri, . . . ,rfc]. 
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Example 2.8. Continuing Example \2.7\ we have that the generators of the radical 
Rad{A) have coordinates 

r., = [1,-3,0,2,0], n = [0,-4,1,3,0], rs = [0,-3,0,2,1] 

in the basis B ~ [1, X2, a;ia;2, x^] . 
We set 

ri = [1, 0, 0, 0, 0], ra = [0, 1, 0, 0, 0]. 



We perform the change of basis to the two multiplication matrices M^^ and 
and obtain: 






1 








■ 




' 3/2 


-1/2 


-3/2 


1 





1 





-1 





1 




-1/2 


3/2 


1/2 














10/3 


-2 


1/3 


and 








-1/3 


1 


-1/3 








5 


-3 


1 










-8/3 


3 


-2/3 








-7/3 


2 


2/3 . 










4/3 


-1 


4/3 



respectively. 

We then have that the multiplication matrices for x\ and xi in yl/i?ad(A) in 
the basis [l,Xi] are 



and — 



3/2 -1/2 
-1/2 3/2 



The eigenvalues of these matrices give the solutions to the system. 



3. Clustered roots 

In this section we consider systems with clustered roots instead of systems with 
root multiphcities. We can think of these systems with clustered roots as being 
obtained from systems with multiplicities via one of the following two ways: 

1. by perturbing the coefhcicnts of the system with multiple roots, 

2. by perturbing the multiple roots to obtain clusters. 

Let f be the system with multiple roots and f be the system with clustered 
roots obtained from f by any of the above methods. Denote by A = C[x]// the 
algebra corresponding to the ideal / generated by the polynomials in f . 



Assumption 3.1. Throughout this paper we make the assumption that the basis B 
for A also forms a basis for A. Note that iff is a well constrained system then for 
"almost all" perturbations f our assumption is satisfied, i.e. the set of perturbed 
systems for which it doesn't hold has measure zero in the space of all systems of 
given degrees. 

If we assume that the basis B for A also forms a basis for A then both the 
multiplication matrices and the matrix of traces are continuous functions of the 
coefficients of the polynomials. Therefore, small perturbations in the coefficients 
of f will result in small changes in the entries of the multiplication matrices and 
the matrix of traces. 
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However, in case 2, when the roots are perturbed, the polynomials corre- 
sponding to the clustered system might end up having coefficients very different 
to those of the original system, even if the radii of the clusters were small. In 
this case, if wo compute the multiplication matrices for the clustered system, the 
entries might not be continuous functions of the perturbation of the roots. They 
not only depend on the magnitude of the perturbation of the roots but also on 
the direction of the perturbation. However, as we shall show in Proposition 13.41 
the matrix of traces is always continuous in the roots. The following examples 
illustrates this phenomenon. 



Example 3.2. We consider three examples of a single cluster of size proportional 
to £ around the origin (0, 0) in consisting of three roots. The first two examples 
demonstrate that the defining equations and the multiplication matrices can have 
different accumulation points as e approaches 0, depending on the direction. The 
third example demonstrate that generally the defining equations and the multipli- 
cation matrices are not continuous at s = 0. 

• First, the roots of the cluster are (0,0), (£,£), {2e,2e). The defining equations 
of these points in C[x, y] are given by x^ — 3ex^ + 2e^x = and y = x, and 
the multiplication matrices in the basis B = {l,a;,a;^} are given by 



= My = 



1 
1 

-2e2 3e 



and lim My = 

e— " 



1 
1 




and the primary ideal defining the multiple root is {x^, x — y) . 

The next example has cluster (0,0), (e, 2£), (fle^Ae). The defining equations 

are x^ — 3ea;^ + 2e^x — and y — 2x. Then is the same as above, but 



M.,. 



2 
2 

-4£2 6£ 



and lim A/„ 



2 
2 




and the primary ideal defining the multiple root is {x^, x — 2y) . 
More generally, the third example has cluster (0,0), {e,ce), (2e, de) for some 
c,d€ R. Then the first defining equation is the same as above, and the second 
equation is y = —^^^x^ + (2c— ^)x which is not continuous in e = 0, unless 
Similarly for the multiplication matrix My. However, the matrix of 



d = 2c. 

traces 



R 



3 

+ de 



ce + de 



4 4 



dV 

dV 
dV 



(with respect to the basis {l,y,y^}) is continuous in e 
limit for every choices of c and d. 



and has the same 



Example 3.3. Continuing with Example \2.8l suppose now that instead of having a 
system with common roots [1, 1] of multiplicity 3 and [—1,2] of multiplicity 2 we 
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have a polynomial system with a cluster of three common roots: 

[[1, 1], [0.9924, 1.0027], [1.0076, 0.9973]] 
around [1, 1] and a cluster of two common roots: 

[[-1,2], [-1.0076,2.0027]] 

around [—1,2]. 

Using the multivariate Vandermonde construction (see for example j32j ). we 
obtained the following multiplication matrices for this system, with respect to the 
same basis as for the system with multiple roots: B = [l,xi,X2,XiX2,x1]. 



The norm of the difference between these matrices and the multiplication ma- 
trices {Ip for the system with multiple roots are very large: 135.41 for the matrices 
of xi and 59.54 for the matrices 0/.T2. Entrywise, the largest absolute value of the 
difference of the entries of the matrices is 55.40 for Xi and 21.70 for X2. 

However, the matrix of traces associated to the system with clusters is 



3.8328 X 10"'^ 


9.9997 X 10"^ 


3.0830 X 10"* 


4.1421 X 10"^ 


3.9951 X 10"^ 




3.7919 X lO""^ 


-2.7338 X lO"'' 


1.2303 X lO"'' 


8.2891 X 10"^ 


1.00004 




3.8527 X 10"" 


-2.7463 X 10"^ 


-1.5183 X 10"** 


1.00000 


3.9969 X 10"^ 




7.73947 


21.69983 


-5.97279 


-16.79084 


-5.67565 




-17.94136 


-54.43008 


13.97610 


41.61207 


17.78328 




3.7831 X 10"'^ 


-2.7103 X 10"^ 


1.00000 


-1.4715 X 10"^ 


4.0017 X 10"^ 


3.8527 X 10"'^ 


-2.7464 X 10"^ 


-1.5183 X 10"** 


1.00000 


3.9969 X 10"^ 


-2.22905 


1.06576 


3.00000 


-7.1053 X 10"^ 


-1.2617 X 10"^ 


-3.23468 


-10.77768 


2.47988 


9.67839 


2.85410 


7.73947 


21.69983 


-5.97279 


-16.79084 


-5.67565 



4.99999 0.99240 7.00269 -1.01796 5.01538 

0.99259 5.01557 -1.01777 7.03349 0.97757 

7.00131 -1.01934 11.00943 -5.04274 7.03192 

-1.01900 7.03226 -5.04240 11.07093 -1.04951 

5.01548 0.97748 7.03339 -1.04838 5.03155 



(3) 



and the 2-norm of the difference between this matrix and the multiplication matrix 
R in (0) for the system with multiple roots is 0.147. 

We have the following result for the entries of the matrix of traces R expressed 
in terms of the roots of the polynomial system. 

Proposition 3.4. The matrix of traces R of the system f(x) = with respect to 
B ~ [61, ... , bn] can be expressed in terms of the common roots {zi, . . . , z„} as 

- n -1 " 

.fc=l 



- jj=i 

where bibj{zk) indicates the evaluation of the polynomial bibj at the point Zfc. 



Proof. Assume that {zi, . . . , z^} are the distinct elements among {zi, . . . , z„} in 
V{I) and let Ui be the multiplicity of z,;. Let Qi be the (unique) primary component 
of / in C[x] whose radical P,; is the ideal of all polynomials vanishing at z^, i = 
1,2, ... ,k. Set Ai = C[x]/(5,;. We have then n.i = dime Ai and / = QinQ2n • • • Qk. 
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Also the ideals Qi are pairwise relatively prime, hence by the Chinese Remainder 
Theorem we have 

A'^ Ai®A2®---®Ak. 

We denote also by Ai the image of Ai in A at this isomorphism. Given any polyno- 
mial g, it is immediate that Ai is an invariant subspace of the multiplication map 
Mg and that the characteristic polynomial oiMg on Ai is {t — g{zi))"'^ . This implies 
that the characteristic polynomial of Mg is HiLi ^ .9(zi))"' = IliLi ^ .9(zi))- 
So the trace of Mg is X^iLi 9i^i)- Therefore 

n 

fe=i 

which proves the lemma. □ 

Note: An alternative proof can be given for Proposition 13.41 using the fact 
that the multiplication matrix Mg is similar to a block diagonal matrix where the 
i-th diagonal block is an rii x rii upper triangular matrix, with diagonal entries 
g{zi), i = l,...,k (cf. [31 Theorem 2]). 

The previous result shows that the entries of the matrix of traces are contin- 
uous functions of the roots, even when the roots coincide. In particular, a system 
with multiple roots and a system with clusters obtained by perturbing the roots 
of a system with multiplicities will have comparable matrices of traces. 



4. Univariate Case 

Before we give our method in full generality we would like to describe our algo- 
rithm in the univariate case. The purpose of this section is to demonstrate the 
simplicity and the accuracy of our technique to compute the approximate square- 
free factorization of a univariate polynomial. As we mentioned in the Introduction, 
our method offers a new alternative to other approximate square-free factorization 
algorithms, such as the one in |25j . 

The following is a description of the steps of our algorithm. Let 

fix) = x'^ + aix'^^^ H h ad-ix + ad e C[x] 

be a given polynomial of degree d with clusters of roots of size at most e. The 
output of our algorithm is a polynomial g{x) G C[x] such that its roots are the 
arithmetic means of the roots in each cluster, with a precision of order of magnitude 

1. Compute the matrix of traces R w.r.t. the basis B = [l,x,x^ , . . . , x"^"^] using 
the Newton-Girard formulas. In this case we have R ~ [si+j]f~ta "^here St 
is the sum of the t-th power of the roots of /. We set sq ~ d and wc find 
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si, ■ ■ ■ , S2d-2 from the coeffrcients of / using the Newton-Girard formulas as 
follows: 

si + ai = 
S2 + aisi + 2a2 = 



Sd + H h aa-isi + da^ = 

Sd+i + aiSd H 1- adSi = 



S2d-2 + 0,lS2d-3 H 1- a-dSd-3 = 0- 

2. Gaussian elimination with complete pivoting (GECP) is used on the matrix 
R until the remaining entries in the partially row reduced matrix Uk are 
smaller than a preset threshold (see Propositions 15.71 and 15. 8p . The number 
of iterations performed, k, is the numerical rank of the matrix R. 

3. Compute a basis of the nuUspace N of the first k rows of the matrix Uk 
obtained after k steps of the GECP. We identify the vectors in N by polyno- 
mials, by combining their coordinates with the corresponding basis elements 
of B. 

4. The smallest degree polynomial in N is the approximate square-free factor 
g{x) of f{x). Its roots are the arithmetic means of the roots in each cluster 
modulo (see Proposition l7.5p . In the case when the matrix R has numerical 
rank d then we take g{x) — f{x) as the square-free factor. 

Example 4.1. (1) Consider the approximate polynomial 

f = {x-{z + 5ie)){x - (z + 52e)){x - [z + S^e)) 
obtained by perturbing the roots of the polynomial 

{x - z)^ = x^ - Zx^z + 'ixz^ - z^. 

Using the basis B = [l,x,a;^] we obtained the matrix of traces R, for which 
the U matrix in the LU factorization obtained by GECP is 

e2^-2<Si53-2a2<S3-2aia2 + 2i5i+2i5| + 2i5|) e^*2 2+^^*2 3 

V 2(-5i53-42«3-'5l«2 + <5i+«2+*3) J 

where are polynomials in the S 's and z 's. 

Using the bound from Proposition \5.8\ for the numerical rank, we have that 
the approximate radical will be defined using the nullspace of the first row of R. 
We obtain the following basis of the approximate radical, 

2 2 2ze [33 + 52 + 5^) + e'^ [Sj + S'i + 51) e {63 + 52 + 5i) . 
{x -z i '-,x-z ^ 
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We choose the element of smallest degree to be the approximate square-free 
factor of f , which is here 

e (S3 +52 + Si) 



We can see that in this case the roots of this polynomial correspond precisely 
to the arithmetic mean of the three clustered roots. 

(2) Consider the approximate polynomial 

f{x) = {x + (-0.98816 + 0.01847/))(2: + (-0.98816 - 0.01847/)) 
(x - 1.02390)(2- - 1.98603)(a; - 2.01375) 

which is a perturbation of the polynomial 

- 7£C* + 19a:^ - 25a;^ + 16a: - 4 = (x - l)^(a; - 2)^. 



The matrix of traces corresponding to f is 



5 

7.00001 

11.00013 
19.00089 
35.00425 



7.00001 
11.00013 
19.00089 
35.00425 
67.01631 



11.00013 
19.00089 
35.00425 
67.01631 
131.05456 



19.00089 
35.00425 
67.01631 
131.05456 
259.16598 



35.00425 
67.01631 
131.05456 
259.16598 
515.47172 



The Uk matrix obtained after 2 steps of GECP on R is 



U2 = 



515.47172 







35.00425 
2.62296 






131.05456 
2.10058 
0.0024342 
0.0029279 
0.0011698 



259.16598 
1.40165 
0.0029279 
0.0035326 
0.0014044 



67.01631 
2.44912 
0.0011698 
0.0014044 
0.00056307 J 



By taking the nullspace of the first two rows of the matrix U2, we obtain the 
following basis of the approximate radical, 

{x^ - 15.01431X + 14.01921, x^ - 7.00397x + 6.00539, 
x^ - 3.00074a; + 2.00102}. 

The approximate square-free factor of f is then 

x^ - 3.00074a; + 2.00102 = {x - 1.00028)(a; - 2.00047). 

We can see that the roots of the output are close to the means of the clusters, and 
the differences are 0.00058 awrf 0.000200 respectively, which are of the order of the 
square of the cluster size (bounded here by 0.03^. 

We refer to the papers of [JH HTJ [51] for other methods that study ap- 
proximate square-free factorization using approximate gcd computation. 



5. LU decomposition of the matrix of traces 

Since the polynomial system with clusters, obtained by perturbing the system with 
multiplicities, has only simple roots, the matrix of traces has full rank. However, 
we can try to find its numerical rank. We will argue below that we can define the 
numerical rank in such a way that it will be equal to the rank of the matrix of 
traces of the corresponding system with multiplicities. 
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In this paper we primarily study the Gaussian chmination with complete piv- 
oting (GECP) [20] in order to estimate the numerical rank and find the numerical 
nullspace of the matrix of traces. However we we will also infer that the singular 
value decomposition (SVD) in our case works similarly to the GECP. 

We would like to note that rounding errors can sometimes result in a matrix 
which is close to a singular one, but where all the pivots are large (see Kahan's 
Example 5 in [23]). This example shows that GECP can be a poor choice for 
numerical rank computations in the presence of rounding errors. On the other 
hand; algorithms for the accurate computations of the SVD of certain structured 
matrices, including Vandermondc matrices, use improved versions of GECP as 
subroutines [151 116| . In our case we prove that the structure of the matrix of traces 
guarantees that we will obtain small pivots which arc proportional to the square 
of the size of the clusters and can therefore use GECP for rank determination. 

We will also show how our results for the GECP of the matrix of traces R 
relate to the singular values of R. In particular we will obtain asymptotic bounds 
for the singular values of the matrix R. Such bounds arc similar to the ones for 
the entries of the Uk matrix obtained after k steps of GECP on R, more precisely, 
we also obtain in this case that the "almost zero" singular values arc proportional 
to the square of the size of the clusters. 

First we study the properties of the Gaussian elimination in the approxi- 
mate setting. We use the following notation for different versions of the Gaussian 
elimination algorithm: 

Definition 5.1. The version of Gaussian elimination in which at the i-th step we 
always select the entry at position {i,i) for pivoting will be referred to as regular. 
We call an m X n matrix M regular if for k := rank(A'/) the first k steps of the 
regular Gaussian elimination on M do not encounter zero pivots. 

Note that GECP on the matrix M computes two permutation matrices P 
and Q of sizes m x m and n x n, respectively, such that for the matrix P M Q the 
regular Gaussian elimination works as GECP. 

In the rest of this section we give results which compare the GECP applied 
to the matrices of traces of the perturbed system and to the system with multiple 
roots. Let Rq be the matrix of traces of the system with multiple roots and let R 
denote the matrix of traces of some perturbation of it. Assume that rank(_Ro) = k. 
Our next result guarantees that for sufRcicntly small clusters, the first k steps of 
the GECP applied to R computes permutation matrices P and Q which make the 
matrix P RqQ regular. 

Proposition 5.2. Let M be an n x n matrix with entries polynomials in x ~ 
[xi, . . . ,xn] over C. Fix z = [zi, . . . , z^] £ , denote Mq :~ A/|x=z> (md as- 
sume that rank(Mo) = k. Then there exists an open neighborhood V of z in 
such that for all points z = [zi, . . . , zj\j] € V if P and Q are the permutation ma- 
trices corresponding to the first k steps of the GECP on the matrix M :— A/|x=i, 
then the matrix P Mq Q is regular. 
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Proof. We call a pair (P, Q) of n by n permutation matrices good if P Mq Q is 
regular, otherwise the pair is called bad. For each bad pair we define an open 
neighborhood Vp.g of z G as follows: For some i < k assume that the regular 
Gaussian elimination on P Mq Q encounters a zero pivot for the first time in the 
i-th step, causing (P, Q) to be a bad pair. Denote by Uq the partially reduced form 
of P Mo Q after the i — 1-th step of the regular Gaussian elimination. Denote by S 
the set of indices (s, t) such that s,t > i and the (s, t) entry of Uo is non-zero, and 
by T the set of indices {s,t) such that s,t > i and the {s,t) entry of Uq is zero. 
Since the rank of P Mq Q is fc, S is non empty. 

Let U be the partially reduced matrix obtained from P M Q via the first 
i — 1 steps of regular Gaussian elimination. Note that the entries of U are rational 
functions of the entries of M and the denominators of these are non zero at z, 
hence are continuous functions of the points [zi, . . . , z^r] in a sufficiently small 
neighborhood of z. In particular, in an open neighborhood U oi z the first i — 1 
steps of regular elimination can be carried out. 

Let the open neighborhood Vp,Q C U C of z be selected such that for 
all [zi, . . . , zn] G ^p.q the entries in T oi U := U\x=z are all strictly smaller 
in absolute value than any of the entries in S of U. By continuity, such open 
neighborhood of z exists, since the required inequalities hold for Uq. 

Finally define V := Vp.g. This is also an open neighborhood of z 

(P,Q) is bad 

since the set of permutations is finite. We claim that for any fixed [zi, . . . , zn] G V, 
if (P, Q) is the pair of permutation matrices corresponding to the first k steps of 
the Gaussian elimination with complete pivoting on the matrix M then (P, Q) is 
a good pair. This is true since PMQ has the property that for each i < k after 
i — 1 steps of the Gauss elimination the (i, i)-th entry of the corresponding matrix 
is maximal in absolute value among the entries indexed by (s, t) ^ (i, i) such that 
s,t > i. But then the (i,i)-th entry in the matrix Uq defined above cannot be 
because of the definition of V. This proves the claim. □ 

In the rest of the paper we will assume that the size of the clusters is a 
parameter e. More precisely, in the following definition we formally explain the 
mathematical setting where our results will hold: 

Definition 5.3. Let z^ — . . . , Zi,m] S C" for i ~ I, . . . ,k, and consider k clus- 
ters Ci, . . . ,Ck of size \ Ci\ = Hi such that^^^^ n.i = n, each of radius proportional 
to the parameter e in the oo-norm around Zi, . . . , z^.' 

C, ={[zi,i + <5,,i,ie,- ■ ■ -l-<5,,i,me], . ■ ■ , 

■ ■ • -l-5,,„,,ie, . . . + 5,, „-,,„£]} (4) 

— {Zj + 5i i€, . . . ,Zi + Si n^e}, 

where \Sij^r\ < 1 for all i — I, . . . , k, j = 1, . . . , n.i, r = 1, . . . , m. Let Uk be the 
partially row reduced form obtained by applying k steps of the GECP to the matrix 
of traces R corresponding to Ci U • • • U Cfe. Then R and Uk have entries from the 
field C{e). 
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Assumption 5.4. Based on Proposition \5.8[ we will assume that if the GECP ap- 
plied to R produces the permutation matrices P and Q then the matrix P RqQ is 
regular, where Rq ~ i?|£=o- To simplify the notation for the rest of the paper we 
will assume that Q = id, i.e. the rows and columns of PRQ = PR correspond to 
the bases 



respectively, where a is the permutation corresponding to the matrix P . This as- 
sumption does not constrain the generality since we may rename B in the definition 



With the assumption that P Rq has rank k and is regular, we can assume 
that aU the denominators appearing in the entries of Uj; are minors of R which are 
non-zero at e = 0. Therefore we can take their Taylor expansion around £ = and 
consider them as elements of the formal power series ring C[[e]]. In this ring we 
shall work with residue classes modulo e^, i.e., in some considerations we factor 
out the ideal (e^) of C[[e]]. 

The results in the rest of the paper are all valid modulo in the formal power 
series setting described above. In practice what this means is that the method we 
propose works up to a precision which is the double of the original size of the 
clusters. 

Remark 5.5. In Dcfinition l5.3l we assume that the clusters arc linear perturbations 
of a set of multiple roots. Note that not all multiplicity structures can be obtained 
as a limit of such clusters with linear perturbation of fixed directions 6i,j as e 
approaches 0. However, as we have seen in Proposition 13.41 the matrix of traces 
at e = is independent of the directions 6i_j, and in fact does not depend on the 
multiplicity structure of the roots. Since all the subsequent results in the paper 
only depend on the matrix of traces and are only valid modulo e^, we do not limit 
the generality by considering only linear perturbations. This is not true however 
for the multiplication matrices, which depend on the multiplicity structure of the 
roots at e = 0, as seen in Example 13.21 

In order to describe the structure of the matrices in the LU decomposition 
of the matrix of traces obtained by GECP in terms of the elements in the clusters, 
we need the following definition: 

Definition 5.6. Let B = [bi, . . . , &„] G C[a;i, . . . ,Xm]", and let zi, . . . ,Zr G C™ be 

not necessary distinct points. We call the n x r matrix 



the Vandermonde matrix o/zi, . . . , w.r.t. B. Note that if r ~ n then the matrix 
of traces in Definition and the Vandermonde matrix are closely related: 



aB = [6^(1) , . . . , 6ct(„)] and B = [bi, . . . ,bn], 



(5) 



ofR. 




R = VV'^. 
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The following proposition gives asymptotic bounds for the entries of the ma- 
trix obtained from a partial Gauss elimination with complete pivoting on the 
matrix of traces R for the case where the n roots of the system correspond to k 
clusters, each of them with rii roots (i = 1, . . . , fc) and radius proportional to e in 
the max-norm. 

Proposition 5.7. Let B = [6i, . . . , 6„] G C[xi, . . . , x™]". Let {zi, . . . , z^} G C™ and 

the clusters Ci, . . . , Cfc around {zi, . . . , Zk} be as in Definition \5.S[ 

Let R be the matrix of traces associated to CiU ■ ■ - UCk and B (see Definition 
\2.4\ and Proposition \ 3.4^ . Let P and Rq := R\e=a be as in Assumption \ 5.4\ and 
assume that P Rq has rank k and is regular. Then, after k steps of the regular 
Gaussian elimination on P R we get a partially row reduced matrix Uk, such that 
its last n — k rows satisfy 

= crrril V-Zt for i = k + l,...,n. (6) 

[cij£ + h.o.t.[£) e C[[eJ] if J > k 

The values of Ci_j G C depends on n, {zi, . . . ,Zfc}, {5s,t\ cind B (we will give a 
bound for Ci^ in Provosition 1 5. 8\] . Here h.o.t.{e) denotes the higher order terms in 
e. Moreover, the formal power series in (0) are convergent in a sufficiently small 
neighborhood of e = 0. 

Proof. To simplify the notation, denote R = P R. The proof is based on the fact 
that after k steps of the regular Gaussian elimination on R, the partially reduced 
Uk matrix has elements (i, j), for i, j = fc + 1, . . . , n, of the form 

det 

-4^ — - (7) 

det(i?W) 

where i?^*^' is the fc x fc principal submatrix of R and r!^^^^ is the (fc+ 1) x (fc + 1) 
submatrix of R corresponding to rows {1, . . . , fc, i} and columns {1, . . . , fc, j}. This 
follows at once from the facts that both the numerator and the denominator of 
(8) stay the same during the row operations performed, and the reduced form of 
i?!^^^' is upper triangular. 

Let V be the n x n Vandermonde matrix of Ci U • • • U Cfe with respect to 
B and recall that R = VV'^, thus R = {PV){V^). Let a be the permutation 
corresponding to P and let aB be as in ([5]). Observe that 



p(fe + l) _ -rrT 

Rij - VaB, Vb 



where V^Bi and Vb are the (fc + 1) x n Vandermonde matrices corresponding 
to Ci U • • • U Cfc and respectively to aBi := [60.(1), ^dr(fe)i ^CT(i)]' ^^"^ ^3 
[61, . . . , 6fe, bj]. Therefore, by the Cauchy-Binet formula we have 

det(i?2+'))= det(KB,,/)det(FB^,/), (8) 

\I\=k+l 

where V^Bij denotes the (fc + 1) x (fc + 1) submatrix of Va-Bi with columns corre- 
sponding to the points in / C Ci U • • • U Cfc, and the summation is taken for all / C 
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CiU- ■ - UCk such that |/| = k+1. Note that aU the determinants in ([8]) are polyno- 
mials in e. Since rank (Fle^o) ~ k, wc have det {VaBi,i)\e=o = dct (VBj,/)|e=o — 0, 
thus they are divisible by e for all i = k + 1, . . .n and / C Ci U • • • U Cfc with 
\I\ = k + 1. Therefore we get that det (J?,^^'*'^^) is divisible by e^. 

Finally we note that the assumption that P Rq has rank k and is regular 
implies that 

det(i?W|,=o)^0, 

which proves that the Taylor expansion of the ratio in (|11[) around e = has 
zero constant and linear terms, as was claimed. The formal power series in ^ are 
convergent in a sufficiently small neighborhood of £ = 0, since they are the Taylor 
series of rational functions with non-zero denominators at e = 0. □ 



From the previous results it follows that if we have k clusters of size Ui , with 



i — 1, . . . , fc, J2i=i 
we get the matrix 



n, then after k steps of GECP on the matrix of traces i?, 



Uk = 



[Uk]i,i 








Ck+l,k+lS 



[Uh]l,n 



Ck + l,n£ 



+ h.o.t.{e) 



(9) 



c„,fe+ie2 

where the constant term in e of [C/fc]i_i is non-zero for i < k. 

The next proposition gives a bound for the coefficient Cij of in It 
also gives an idea of the magnitude of the threshold one can use to decide on the 
numerical rank which would additionally indicate how small the size of the clusters 
need to be for our method to work. 



Proposition 5.8. Let B = [&i,...,6„] £ 
Let the clusters Ci, . . . , Ck around {zi, . . 



* = 1,.. 

to Ci U 



,k, j — 1, . . . , n,; 
• U Ck and B. 



1,. 



C[xi,...,XraT. Let {zi,...,Zfc} G C". 

. jZfc} he as in ^ with \5i_j^r\ < 1 for all 
Let R he the matrix of traces associated 



Let b' he such that 



— (z ) 

dXr 



(10) 



Assume that the GECP applied to R also implies complete pivoting on R\^^q. 
Then the hound for the coefficients Cij of in the Uk matrix, ohtained after k 
steps of the GECP applied to the matrix of traces R, is given hy 

k,i <«-(br, 

where a = 4(n — k){k + 1)^to^. 
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Proof. We denote R := PR, where P is a permutation matrix such that the 
first k steps of GECP applied to both PR and P R\e=o is well defined and the 
same as regular Gaussian elimination. Note that we need the assumption that 
GECP applied to R also implies complete pivoting on R\^^q since Proposition 
15.21 only implies that P R\e=o is regular, but below we will also need the pivots 
in P R\e=o to have maximal absolute values. One can achieve this by making the 
right selection among equal possible pivots while performing GECP on R. We will 
use this assumption at the end of the proof. 

Denote the bases corresponding to the rows and columns of R by aB = 
[^<t(i), ■ • ■ , ^<t(«)] and B = [bi,..., 6„] as in (O. 

The partially reduced Uk matrix has elements for i, j = fc + 1, . . . , 71, of 

the form 

det 

— (11) 

det(i?W) 

where R^'^^ is the fc-th principal submatrix of R and Rg^gn is the (fc + 1) x (fc + 1) 

submatrix of R corresponding to rows B' := [fccr(i), • ■ • , &cr(fc) i &cr(i)] and columns 
B" := [61, . . . , bj]. In order to get an upper bound for \cij\, we will get an upper 

bound for the coefficient of in det {Rg^g,, ) and divide it by the constant term 

of I det (^e^'))]. 

Fix i, j e {1, . . . , n}. We will use the Cauchy-Binet formula 
det(i?^'+il) = det{VB',i)det{VB",i), 

\I\=k+l 

where the summation is for / C Ci U • • • U of cardinality fc + 1, and Vb'j 
and Vb" J are the Vandermonde matrices corresponding to / w.r.t. B' and B", 
respectively. Since the derivative of the determinant of a matrix is the sum of 
determinants obtained by replacing one by one the columns by the derivative 
of that column, after expanding the determinants in the sum by their columns 
containing the derivatives, we get 



a det (Vb'j) 



de 



db' 



E=0 



^ ^± (II'^M^(z)|e=0 I det(l/B,_{b,}j_{,})|_Q, 



b'eB' zei 



where t is the coefficient of e in the <-th coordinate of z e /. We can obtain a 



similar expression for 



adct(yB„.,) 



£=0 



Note that det ( VB/_{b},/-{z})|,_o non-zero only if I\e=o = {zi, • . • ,za;} U 
{zi} for some i = 1, . . . , fc and zj^^o = z.^. In that case / — {zjlg^o ~ {zi, . • . , z^.}, 
which we denote by 

3 := {zi, . . . ,Zfc} 

for simplicity. 

Thus we have that if /|e=o ~ {zi, . • . , z^.} U {z.^} then 
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ddet{VB',i) 



de 



< E E |det(^B'-m,3)lE 



db' 
dxt 



(z) 



6 = 



using that \5^ t \ < 1- Therefore, we get that the coefficient of in det [R^i g,,) is 
bounded by 

4 |det(Fs'-{b'},3)det(FB,,_{6"},3)| 
/ 



b" £B'' 



E E 

1 l/l=fe+l \t=l 



(z.) 



E 



96" 

9a; t 



(^0 



using the fact that there are two possible ways to pick z G / with zj^^o = from 

/|e=0 = {zi, . . . ,Zfc} U {zj. 

Using the upper bound b' and counting the number of times we can choose 
/ C Ci U • • • U Cfe such that |/| = A; + 1 and /|e=o = {zi, ■ • ■ , z^} U {z^}, we get 



Am^h' {n-k) ^ 



^ |det {VB'-{b'},?,) det (l^i3"-{b"},3) | 



On the other hand, wc have that if RB'-{b'},B"-{b"} is the matrix of traces 
with rows corresponding to the B' — {fe'} and columns corresponding to B" — {6"} 
then 

Acl{RB'-{b'},B"-{b"}) = ( ) dot (yB'-{b'},3) dot (yB"-{6"},3)- 



Therefore, the bound for the coefficient of in det {Rb/^b") 



Am^h' (n-k) ^ dot (^B,_{b,}.B/,_{f,/,}) 



b' eB' 
b" es" 



£-0 



(12) 



Next we use the assumption above on i?|e=o = PR\e=o to have maximal 
pivots in the first k diagonal entries to get 



det (i?('^')|e=0 > I det 



6 = 



(13) 



which is true since the left and right hand side of divided by | det {w-''~^^)\e=o 
give the absolute values of the entries of the partially row reduced matrix after 
fc — 1 steps of GECP. Therefore we can replace | det (_RB/_{f,/} b"-{6"})|6=o by 
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I det (^('=))|e=o in (US) and divide the expression ^ by | det (^('=')|e=o, thus we 
get the following bound for the coefficient Ci,j of e"^ in the Uk matrix: 

□ 

Remark 5.9. The above proposition gives estimates in terms of {zi, . . . , Zfe}. which 
we do not assume to know a priori. The following heuristic methods can be used to 
check whether the estimated numerical rank is correct, given a required precision 
e. Assuming that we know the magnitude of the coordinates of the roots, we may 
compute the matrix of traces corresponding to n randomly chosen distinct roots 
which have the same order of magnitude as the original roots. Then comparing 
the diagonal entries of the {/-matrices obtained by applying the GECP for the 
matrices of traces, we can set the numerical rank to be the first entry where 
the discrepancy is of order e"^. Another heuristics is to increase k one by one, 
compute the approximate radical ideal (see Definition 17. ip corresponding to the 
case when the numerical rank of R is k. Compute the roots of the approximate 
radical ideal, and substitute them back into the original system. If the error is 
of order of magnitude e, accept k and the computed approximate radical as the 
output. 

Example 5.10. Continuing Example \3.3\ we apply the GECP to the matrix R 
defined in (0). After two steps of GECP we obtain the following matrix: 



U2 = 



11.07093 -5.04240 7.03226 -1.01900 -1.04951 

8.71265 2.18381 6.53716 6.55387 

0.454213 X 10"* 0.7407 x lO"'"* 0.178036 x lO"''' 

0.7397 X 10"^ 0.728 x 10"'' 0.41955 x 10"* 

0.188071 X 10"^ 0.52002 x 10"* 0.657084 x 10"'' 



with columns permuted so that they correspond to the basis [a;ia;2, 2:2, Xi, 1, a;^]. 
Note that the largest entry in the 3x3 bottom right corner of U2 is between e and 
(here e w 0.01 in this example). Thus we consider the numerical rank of R to 
be 2. From the nullspace of the first two rows of U2 we can obtain the following 
approximate multiplication matrices: 

, _ r 1 1 A^' - r 1.49973 -0.49972 

"=1 ~ [ 1.00382 -0.37849 x 10"-'' J "2 ~ [ -0.5016325 1.50162 

(see Section^ below for more details on approximate multiplication matrices). 
The eigenvalues of A4'^ and M'^ are respectively 

{1.000018,-1.003803} and {0.9999943,2.001349}. 

Note that these eigenvalues are close to the avarages of the coordinates of the 
roots in the two clusters. 



6. Singular Values of R 

Using the previous results we will now study the singular values of the matrix of 
traces i? of a system with clusters of roots. We denote R P R, where P is a 
permutation matrix obtained by k steps of GECP applied to R and we assume 
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that P R\e=o is regular, as in Assumption 15.41 Let Uk be the matrix obtained after 
k steps of GECP on the matrix R, as in ([9]). Let Uk he the matrix obtained after 
replacing the last n — k rows of Uk by zeros. Let be such that R ~ L^Uk (in 
other words is the transformation matrix obtained after k steps of GECP on 
R). Let R = LkUk. Using the submultiplicative property of matrix norms, we have 
that 

\\R-R\\f = \\LkUk - LkUkllF < \\Lk\\F\\Uk - Uk\\F 
where j| • \\f denotes the Frobenius matrix norm. 

Let ci > ■ • • > cr,i be the singular values of R, which are also the singular 
values of R. Since by definition cr,; is the 2-norm distance from R to the nearest 
rank i matrix, and R is an n by n matrix of rank fc, we have that 

CTn < ■■■(Jk+l < ||^-i?||2- 

Given that the 2-norm of a matrix is smaller than or equal to its Frobenius 
norm, we have 

an < •••CTfe+l < \\Lk\\F\\Uk-Uk\\F. 

Since we are using GECP it is easy to see that 

[ifc],j < 1 for £dlj = l,...,k,i>j 
and the matrix obtained after k steps of GECP is of the form 



1 











Jfe+i,fc 



Therefore we have 



\L 



k\ F 



< 



2n + 2nk - k^ - k 



From Proposition 15.81 we have that for i,j = k 
are of the form 

[Ukl,j ^ije^ + h.o.t.{e), 

where uj — 4(n — k)(k + l)^m^b' and b' is defined in ([10 
We therefore have 



1 . . . n, the elements of Uk 



\\Uk-Uk\\ 



k\\F 



\ 



E (|tM.j)^<(' 

ij=fe+l 



k)ujs'^ + h.o.t.{s). 
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We summarize the above argument in the next Proposition, showing that the 
k + 1-th singular value of R is asymptotically equal to e^. 

Proposition 6.1. Let R be the matrix of traces associated to Ci U • • • U and 
B where the clusters Ci, . . . , C^. around Zi, . . . , G C™ are as in ^ and B = 
. . . ,6„] £ C[a:i, . . . jX^]". Let (Ti > • • • > (t„ 6e the singular values of R. Then 

cTfc+i = Q.e'^ + h.o.t.{e) 

where 



2n + 2nk — k'^ — k , 2 



n <4(n-fc)^(fc + l)^m^y ■ (b')^ 

and b' is defined in ilO\) . 

Example 6.2. Continuing Examvle \3.3\. we compute the singular values of the ma- 
trix R defined in 

[22.8837, 14.2433, 0.448334 X 10"^ 0.174904 x lO"'^, 0.594796 x 10"^]. 

We have that the third singular value is between e and (in this example e ~ 0.01 ), 
thus we can set the numerical rank of the matrix R to be 2. Note that the 2-norm 
distance of the matrix R from i?|E=o ^s not the same order of magnitude as the 
third singular value, it is 0.147 as was computed in Example lS.Sl This is the reason 
why we used the partial LU-decomposition of R and not R\e=o to obtain a bound 
for ak+i- 



7. Approximate Radical Ideal 



Using our previous results, we can now define the concept of an approximate 
radical ideal and describe its roots in terms of the elements of the clusters. 

Definition 7.1. Let B ^ [fei, . . . , 6„] € C[xi, . . . , Xm]"' and the clusters Ci, . . . ,Ck 
be as in Definition \5.3[ Let R be the matrix of traces associated to Ci U • • • U Cfc 
and B. Let the permutation matrix P corresponding to the permutation a obtained 
after k steps of GECP on R as in Assumption \5.4\ so that the rows and columns 
of R :~ P R correspond to uB and B, respectively, as in We define the vectors 
V,: 



ii^j G C(e)'^ for i = 1, . . .m and j = 1, . . . k, as the solutions of the following mk 
linear systems: 

i?(^V,j=r,j z = l,...,m j = l,...,k, (14) 

where the left hand sides are always the k x k principal submatrix of R, while for 
any fixed i and j the right hand side of ([I^ is defined as 

Tr{x^bjb„(i)) 



Tr{xibjb„(k)) 



(15) 



Note that one can compute the vectors r^ j the same way as the columns of the 
matrix of traces. Then we define the following mk polynomials: 
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:= - |^^[v,j]s6sj i = l,...,m, j = l,...,k. (16) 

We will call the approximate radical ideal of the clusters Ci U • • • U Cfc the 
ideal generated by 

^ ■= (/ij ■ i = l,---,m, j = l,...,k). 

We also define the approximate multiplication matrices of the radical o/ Ci U • • • U 
Cfc with respect to the basis [hi, ... , bk] to be the matrices M^^, . . . , M^^^ E C(e)'^^'^ 
where 

[M'^Jj,s ■■= Wuj]s i = l,...,m, j, s 1, . . . , fc. 

Remark 7.2. We can also define the approximate multiplication matrices of the 
radical of Ci U • • • U from a system of multiplication matrices of Ci U • • • U Cfe 
with respect to B by changing the basis as follows: Let rk+i, ■•■,?'„ G C(e)" be 
a basis for the nuUspace of the first k rows of P R. Let wi , . . . , Wfe € C" be such 
that B' := [wi , . . . , Wfc , rk+i , . . . , r„] forms a basis for C(e)" . Let Mx^ , ■ ■ ■ , G 
C(£)"^" be the multiplication matrices of the clusters Ci U • • • U C/t with respect 
to the basis B' . Then the approximate multiplication matrices of the radical of 
CiU- • -UCfe with respect to [ui, . . . ,Vk] arc the matrices M^^, . . . , e C(£)''^'' 
obtained as the principal fc x fc submatrices of M^^ , • ■ • , Me„ , respectively. Note 
that the eigenvalues of M^^ are the Xi coordinates of the elements of the clusters 
reordered in a way that the first k correspond to one eigenvalue from each clus- 
ter. However, we also remark that we have to be careful with the multiplication 
matrices Af^^^, . . . , Mx^ since they are not always continuous at e = 0, as noted 
in Rcmark l5.5[ thus we cannot consider their entries as elements of C[[£]]. That is 
the reason we chose to define the approximate radical as in Definition 17. II 

The next proposition asserts that when e — Q our definition gives the multi- 
plication matrices of the radical ideal. 

Proposition 7.3. Using the assumptions of Definition \7.1\ the coordinates of the 
vectors vij E C{e)'' defined in are continuous in e = for all i = 1, . . . ,m 
and j ~ l,...,fc. Furthermore, the points Zi,...,Zfc are common roots of the 
polynomials {/ij \s=o ■ i = 1, . . . , to, j = 1, . . . , k}, and the matrices 

form a system of multiplication matrices for the alge bra C[x]/V7. 

Proof. Using Assumption 15.41 the continuity of the coordinates of the vectors 
e C(e)'^ follows from our assumption that the fc x fc principal submatrix W-^^ 
of R is nonsingular at e = 0. 
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Next we show that zi , . . . , are roots of /; j |e=o for all z G {1, . . . , m} and 
j e {I, . . . ,k}. Fix i and j. Assume that 

Xibj - ^^Wij^s^s^ =0 (17) 

is satisfied by zi, . . . , z/;, which is equivalent to the column vectors 

Wjj := [wijs, ■ ■■ , Wij^k, -1]^ (18) 

satisfying the homogeneous linear system with coeSicient matrix W defined to be 
the transpose of the (k + 1) x k Vandcrmonde matrix of Zi, . . . , with respect to 
[bi, . . . ,bk,Xibj]. 

On the other hand, by (|14p . the vector [v^j j — l]e=o is in the nuUspace of the 
A; X (fc + 1) matrix [R^'''>\r.ij]^=Q. We have 

where Vi and V2 are the Vandermonde matrices of Ci , . . . , Cfc at e = correspond- 
ing respectively to [&cr(i), ■ • ■ , ^'cr(fe)] and [bi, . . . ,bk, Xibj], thus is the same as 
W except the row corresponding to Zs is repeated Us times for s = 1, . . . , A:. This 
implies that the nuUspace of is a subset of the nuUspace of [w''^rij]^=Q. But 
since both nuUspaces has dimension one, we must have w^j = | — l]£=o- i-C- 
fi,j\e=o = is satisfied by zi, . . . ,Zfc. 

Next we prove that the matrices M^^^ \e=Oi ■ • ■ i \^^q form a system of mul- 
tiplication matrices for C[x]/a/7. First note that for any 5 G C[x], if z is a common 
root of the system 

gbj - ^ Cj^sbs = j = l,...,k 

s=l 

and z is not a common root of 5i, . . . , 6^ then g{z) is an eigenvalue of the matrix 
■= [cj,s]|,s=i ^ith corresponding eigenvector [bi{z), . . . ,bk{z)]'^ 7^ 0. Our as- 
sumption that = has rank k implies that the vectors [61 (z^), . . . , 6fc(zs)]^ 
for s = 1, . . . ,k are linearly independent, thus they form a common eigensystem 
for the matrices M^Jg^oj ■ • ■ i le=o- Thus, they pairwise commute and their 
eigenvalues are the coordinates of Zi, . . . , z^, proving the claim. 

□ 

Remark 7.4. Without further assumptions on the polynomials 61 , . . . , 6fc we cannot 
guarantee that the polynomials fij\e=o have no roots outside of zi, . . . ,Zfc. For 
example, if fc = d = 1 and Zi = c 7^ but &i = x, then /n — — cx which 
also have as a root. However, if we assume that &i, . . . , 6^ have no common roots 
in C™ (e.g. 1 € {&!, . . . ,bk}) then all common roots of the polynomials fij\e=o 
correspond to eigenvalues and eigenvectors of M^.\^=o- Since Zi,...,Zk already 
provides a full system of eigenvectors for M^Je=Oi the polynomials fij\e=o cannot 
have any other distinct root. 
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Our last result gives an asymptotic description of the roots of the polyno- 
mials {fij} in the case when e ^ 0. Since the coordinates of the vectors v^j- are 
continuous in e = we can take their Taylor expansion around e = and consider 
them as elements of the formal series ring C[[e]]; as described in Definition 15.31 
In this setting we will show that the roots of the system {fij} are the centers 
of gravity (or arithmetic means) of the clusters, modulo e^. Since the arithmetic 
mean of a cluster is known to be better conditioned than the individual roots in 
the clusters (c.f. [30l[T0j), our result is therefore stable for small enough values of 
e. 

Proposition 7.5. Let B ^ [bi, . . . ,bn], {zi, . . . , Zfe} and for i = 1, . . . , fc 

Ci ={lz,,i +<5,,i,ie, . . . + (5i,i.,„e], . . . , 

■ ■ ■ + Si_„.^ie, . . . .Zi,™ + <5,,„^,„e]} 

be as in Definition \5.3[ Let ~ [Cs.i: ■ • ■ i ?s.m] for s = 1, . . . fc he defined as 



-e i = 1, 



(19) 



Then ^i, . . . , satisfy modulo the defining equations {fi.j} of the approximate 
radical ideal of CiU ■ ■ ■ U Ck defined in Definition \7.1\ 

Proof. Fix i e {1, . . . , m} and j ^ {1, . . . ,k}. Define W to be the transpose of the 
(fc + 1) X fc Vandermonde matrix of ^i, . . . , with respect to [bi, . . . , bk, Xibj], i.e. 

W:= [bt{l)\{x,b^)iu'^ 
Also define S to be the (fc + 1) x fc augmented matrix 



S := 



where i?'^'^-' and r^.j was defined in Definition 17. II Assume that 

/ fe \ 



''^^Wij.sbs ] — mod 



(20) 



is satisfied by = [S.s,i,- 
column vector 



■ ,Cs,m] for s = l,...,fc, which is equivalent for the 



Wij := [wij^i,...,Wij^k,-lf (21) 

to satisfy the homogeneous linear system with coefficient matrix W modulo . On 
the other hand, from the definition of the approximate radical ideal in Definition 
17.11 we also have that the augmented vector [vij | — 1] is a solution of the homo- 
geneous system corresponding to S. By our assumption that det(i?^'^^ |e=o) 7^ 0, 
we also have that det(i?'''"-') ^ mod e^, which implies that both S and W have 
nuUspace of dimension 1 modulo e^. Thus it is enough to show that Wj j is in the 
nuUspace of S modulo e^, that will imply that Wi j = [v.; j| — 1] mod e^. 
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Write 

= w'^^ + w^^ e W = W^°'> + W^^h S = S^°'>+S^^h mod 

At e = we showed in the proof of Proposition 17.31 that if wf^ is in the nuUspace 
of then it is also in the nuUspace of 5'°-'. 

It remains to prove that W'^^^v^l"] +W^"^vf[^^ = implies S^^^v^["] +S'^^^vf[^j = 
0. We use the fact that 

and S^''^ = V,W^'^ + {w['^Y V^^ 

where Vi and V2 are the Vandermonde matrices of CiU - • -UCfe at e = correspond- 
ing respectively to [&cr(i), • ■ • , 6cr(fc)] and [bi, . . . , bk, Xibj], Wi is the Vandermonde 
matrix corresponding to ^1, . . . with respect to (&o-(i), ■ ■ • , and Vi and V2 

are the same as Vi and V2 , except the row corresponding to appears only once 
and it is multiplied by Ug . Since wf'j is in the nuUspace of W^'^^ , it is also in the 
nuUspace of , thus it remains to prove that 

V.W^'^wf;^ + V.VM'I = 0. (22) 

Since W^^^w\^j — —W^^^^w^^j by assumption. (p2)) is equivalent to 

-ViW^"^ + ViV^] Mvllj = 0. 

But it is easy to see that ViW^^"'' — V1V2 , which proves the claim. 

□ 

As a corollary of the previous proposition we get that modulo the approx- 
imate multiplication matrices M'^^ , ■ ■ ■ , M'^^ form a pairwise commuting system of 
multiplication matrices for the roots ^1 , . . . , ^a: . 

Corollary 7.6. Using the notation of Definition and Proposition \ 7.5\ we have 
that for all i ~ I, . . . , k and j = 1, . . . , d 

where 



k 



s=l 



Thus the vectors {e^ }^^i form a common eigensystem for the approximate multi- 
plication matrices M^^, . . . , M'^^ modulo , which also implies that they are pair- 
wise commuting modulo , i.e. the entries of the commutators M'^^M'^. — M'^.M'^^ 
are all divisible by . 

Remark 7.7. In practice, for any particular choice of e € 1R+ the system {/ij} is 
not necessary consistent. Also, the approximate multiplication matrices M^^ , ■ . ■ , -M^^ 
are not pairwise commuting, and therefore not simultaneously diagonalizable. 
However, one can take any consistent subsystem of {/i.j} such that it defines each 
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of the coordinates and solve this subsystem in order to obtain the solutions. An- 
other approach is the one described in [301 113 • If the distance of the clusters from 
each other were order of magnitude larger than the size of the clusters then a ran- 
dom linear combination of the matrices M^^ , • ■ • , will have all its eigenvalues 
distinct with high probablility. Using the eigensystcm of this random combination 
one can approximately diagonalize all of the approximate multiplication matrices 
M^^, . . . , M'^^. Then by CoroUarv 17.61 and [10} Proposition 8] the entries outside of 
the diagonal of the resulting matrices will be small, asymptotically e^. Taking the 
i-ih. diagonal entry of these nearly diagonal matrices will give the coordinates of 
the i-th root of the approximate radical, which by Proposition l7.5l is approximately 
the arithmetic mean of a cluster. 



Example 7.8. Our last example is similar to Example \3.3\ but here we increased 
the size of the clusters. Consider the polynomial system given by 

fi =xl+ S.gggSOxiSa - 5.89970ii + 3.81765a;i - 11.25296a;2 
+ 8.33521 

f2=xl + 12.68721a;?X2 - 2.36353j:? + 81.54846x1X2 - 177.31082x1X2 

+ 73.43867x1 - Xa + 6x2 + X2 + 5 
fs=xl+ 8.04041x1X2 - 2.16167X? + 48.83937xiX2 - 106.72022xiX2 

+ 44.00210x1 - xi + 4x2 + X2 -f 3 

which has a cluster of three common roots, [0.8999, 1], [1, 1], [1, 0.8999] around 
[1,1] and a cluster of two common roots, [—1, 2], [—1.0999, 2] around [—1,2]. The 
clusters has size at most e = 0.1. Using Chardin's subresultant method, we ob- 
tained the multiplication matrices for this system, with respect to the basis B = 
[1, xi, a;2, 2;i2;2, a;^] and computed the matrix of traces associated to the system, 
which is 



5 0.79999 6.89990 -1.40000 5.01960 

0.79999 5.01960 -1.40000 7.12928 0.39812 

6.89990 -1.40000 10.80982 -5.68988 7.12928 

-1.40000 7.12928 -5.68988 11.45876 -2.03262 

5.01960 0.39812 7.12928 -2.03262 5.11937 



(23) 



C/2 ■ 



After 2 steps of GECP on the matrix of traces we find the partially reduced matrix 

11.45876 -5.68988 7.12928 -1.40000 -2.03262 

7.98449 2.14006 6.20472 6.11998 

0.01039 0.00799 0.02243 

0.00799 0.00728 0.01544 

0.02243 0.01544 0.06796 

with columns permuted to correspond to the basis [2:1X2, 2:2, xi, l,x^]. 
We also computed the singular values of R: 

[24.06746, 13.29215, 0.04397, 0.00362, 0.00035]. 

We indeed have that the entries in the last three rows of U2 and the third singular 
value (73 are of the order of , which would determine the numerical rank of R to 
be 2. 
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By considering its last three rows of U2 as zero, we compute the nullspace of 
the resulting matrix, which gives the following generators of \fl jl , 
r-A^X2 - 1.46302 + 0.510803x1, 
r-4 = X1X2 + 0.51920 - 1.505323x1, 
rs = X? - 1.01587 + 0.08562x1. 

From these we can define the multiplication matrices for xi and X2 in C[x]/\/7 
in the basis [1, xi] : 

, _ r 1 1 Ayf' - r 1-46302 -0.51080 

=^1 ~ [ 1.01587 -0.08562 J ^2 ~ [ -0.51920 1.50533 

These matrices do not commute but their commutator have small entries: 

, , _ ^1 ^1 _ \ -0.000293 -0.00143 
^rri^xa ^X2^xi " |^ 0.00147 0.000293 

Thus the multiplication matrices are "almost" simultaneously diagonalizable. Fol- 
lowing the method in jlOj . we get the following approximate diagonalizations of 
M'^^ and M'.^^ using the eigenspace of M.'^_^ 

, r -1.05162 0.001765 1 , \ 1.99959 -0.001768 

''I ~ |_ 0.00116 0.966001 J ~ |_ -0.001169 0.968759 

The corresponding diagonal entries give the solutions [—1.05162, 1.99959] and 
[0.966001,0.968759] which are within 0.00167 distance from the centers of gravity 
of the clusters in the oo-norm. 
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